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Abstract. The computation of the symmetric regularized-long-wave (SRLW) equation, 
which describes weekly nonlinear ion acoustic and space-charge waves, is dealt with in this 
paper. The numerical scheme to be proposed applies the Fourier pseudospectral discretiza- 
tion to spatial derivatives in time space, with time advance accomplished in phase space by 
a integrator based on trigonometric polynomials which is fully explicit. Extensive numer- 
ical tests are reported, which are geared towards understanding the accuracy properties 
and studying the collisions of solitary waves in the SRLW. The results suggest that the 
scheme is of spectral-order accuracy in space and second-order accuracy in time. Also, 
some intriguing phenomena in the solitary wave collisions are observed in simulations. 



1. Introduction 
Consider the symmetric regularized-long-wave (SRLW) equation, 

(1.1) utt-u xx + ^(u 2 ) xt -u xxt t = 0, t > 0, 

which was first derived by Seyler and Fenstermacher [13. from a weakly nonlinear analysis 
of the cold-electron plasma equations appropriate for a strongly magnetized nonrelativistic 
electron beam such that the fluid motion is constrained to one direction. This equation was 
shown to describe weekly nonlinear ion acoustic and space-charge waves, and the real- valued 
u(x, t) corresponds to the dimensionless fluid velocity with a decay condition liro.| a .|_ >00 u = 0. 

The SRLW equation is explicitly symmetry in the x and t derivatives, and is very similar 
to the regularized-long-wave equation which describes shallow water waves and plasma drift 
waves [HE]- The existence and uniqueness of solutions for the SRLW have been established, 
e.g., in [H[5]. R is well-known that the SRLW equation (II. ip possesses the solitary waves of 
the form [5|fl3] 

(1.2) u c (x, t) = 3(0 ~ ^ sech 2 ( J ^-^{x -ct)) . 



Note that the two branches of solitary waves for the velocity c in the ranges c > 1 and 
c < — 1 corresponds to the right and the left traveling solitary waves of the same type, and 
hence the SRLW equation describes bidirectional propagation. 
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The SRLW equation (jl.ip has an equivalent form as a first-order system: 

(1-3) u t + p x + - (u 2 ) x - u xxt = 0, 

(1.4) p t + Ux = 0, iEK, t>0, 

from which one can determine four invariants of the motion for the SRLW [5lll3j: 



(1.5) I(u) = I u(x,t)dx, I(p) = / p(x,t)dx, 

J —oo J —oo 

/oo 
(u 2 (x, t) + ul(x, t) + p 2 (x, t)) dx, 
-oo 



oo 



(1.7) V(u, p) = J \u(x, t)p(x, t) + ^u 3 (x, t)\ dx. 

As pointed out in [13], the system (jl.3)) - (jl.4[) dose not possess the Painleve property which is 
a sufficient condition for integrability and (|1.3p - (|1.4j) has only a finite number of polynomial 
conservation laws, all of which are represented by (|1.5p -( fL~7|) . 

In this paper, we are going to address the numerics of the SRLW. Due to its wave- 
type shape and the decay condition, the spectral-type approach has been extensively used 
for the SRLW. For instance, in [HUB] the Fourier spectral/pseudospectral discretization 
was applied to solving the equivalent first-order system (jl.3p - (|1.4j) with the time advance 
accomplished by backwards Euler method. Another category of the numerical methods 
available in literatures for the SRLW is the time-domain finite difference (TDFD) methods. 
Among of them three conservative schemes were proposed and analyzed in [15], which 
conserve some discrete quantities analogous to (jl.5p and (jl.6p . The key point in designing 
conservative TDFD schemes for wave-type differential equations is to use an appropriate 
temporal discretization, like the classical Crank-Nicolson method, and the rigorous proof of 
conservation proceeds by means of 'discrete energy method' with the help of summation by 
parts formula. 

We would like to comment that one can readily obtain the spectral-type conservative 
schemes for the SRLW by combing the Fourier spectral/pseudospectral discretization with 
the finite difference temporal discretization used in conservative TDFD methods. The 
conservative quantities can be proven by analogous arguments as TDFD noting that the 
integration by parts holds in the finite-dimensional functional space spanned by Fourier 
bases and summation by parts also holds in the vector space spanned by discrete Fourier 
bases [7]l9l ll2irT4"] . Although conservative schemes have nice mathematical properties, in 
practice at each step certain fully nonlinear problem has to be solved very accurately which 
requires many iterations in implementation and is quite time-consuming. If the nonlinear 
problem is not solved very accurately, the conservation dose not hold in computation (see 
e.g., [3]). In the view of these observations, the spectral-type conservative scheme for the 
SRLW is not examined in this paper. 

The aim of this work is to propose an efficient and accurate numerical discretization for the 
SRLW equation (jl.lj) . and apply the scheme to study the collisions of solitary waves in vari- 
ous setups. The scheme to be proposed is based on the application of Fourier pseudospectral 
discretization to spatial derivatives, followed by a numerical integration via trigonometric 
polynomials to the resulting second-order nonlinear ODEs for the Fourier coefficients in 
phase space. The numerical integration first writes out the exact solutions of the second- 
order nonlinear ODEs in integral form, and then approximates the integrals by proper 
numerical quadratures such that the resulting scheme is fully explicit. Such integrator is 
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somewhat similar to a Gautschi-type exponential integrator [8l[T01[TT] and a Deuflhard-type 
integrator [6j[TT] , which have been successfully used in the semi-discretization of wave- type 
equations (see, e.g., OHD]) and approximates much better than frequently used finite dif- 
ference integrators under the same time step (see, e.g., [3])- But it also differs from the 
standard Gautschi-type or Deuflhard-type integrator, since the integrands arising in the 
mixed derivative term in (jl.ip suggest a four-level scheme while the standard Gautschi-type 
or Deuflhard-type integrator for wave-type equations results in a three-level scheme. 

The rest is organized as follows. In Section[2]the efficient numerical method is proposed in 
very details. In Section [3] we first test the accuracy properties of the proposed method, and 
then apply it to study the collisions of solitary waves in the SRLW. Finally, some concluding 
remarks are drawn in Section [H 

2. Numerical method 

In this section we shall propose a numerical discretization for the SRLW (jl.ip . In practice, 
we truncate the whole space problem (jl.ip into an interval $7 = (a, b) with periodic boundary 
conditions, i.e., consider the following initial-boundary- value problem: 

(2.1) u a - u xx + - (u 2 ) xt - u xxtt = 0, i>0, 

(2.2) u(x,Q) = u {0) (x), u t (x,0) = u {1) (x), x £Tl := [a,b], 

(2.3) u(a,t) = u(b,t), u x (a,t) = u x (b,t), t>0. 

Such boundary conditions are inspired by the decay condition limM_ >00 u = 0. Choose mesh 
size h := (b — a)/M with M an even positive integer, time step r > 0, and denote the 

mesh points with coordinates (xj,t n ) = (a + jh, nr) for j = 0, 1, . . . , M and n = 0, 1, 

In what follows we will approximate the spatial derivatives in time space along the Fourier 
pseudospectral approximation path, followed by advancing the time in phase space by a 
numerical integration based on trigonometric polynomials. 
Assume 

M/2-1 

(2.4) u(x,t)pt ui(t)exp (im(x - a)) , 

l=-M/2 

for x £ O, and t > 0, with ui(t) the discrete Fourier transform coefficient of the Ith. mode. 

^ = Tl 2 u fo>*w(-^-) - w = &37' l = ~M/2,...,M/2-l. 
Approximating the spatial derivatives in (|2.1|) by Fourier pseudospectral discretization (cf. 

[TKHKEKH]), i.e., 

M/2-1 M/2-1 

-d xx u(x, t) w ^2 ViUi(t)exp (ifii(x - a)) , d x u(x,t) « ^ i/^(t)exp (ifn(x - a)) , 

l=-M/2 l=-M/2 

and noting orthogonality of the Fourier functions, we obtain the following second-order 
ODEs in phase space, for I = -M/2, . . . , M/2 - 1, 

(2.5) iL«,(*) + tfui(t) + y + tf^l(t) = 0- 
Now we discuss a numerical integration for (|2.5p in very details. 
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2.1. Integration based on trigonometric polynomials. For t n (n = 0, 1, . . .) any given 
time and s S R, we rewrite the ODEs (|2.5p as 

d 2 i 

Mtn + S)+ Pfu^tn + S)+ . m 2 J P(S) = 0, 



ds 2 



where, 



A 



1 + /*? 



2(1 + /if)' 



/TOO = («'),(*» + «). 



The analytical solutions of the above second-order ODEs can be written explicitly thanks 

to variation-of-constants formula, 

(2.6) 

rs d 



\~ sin(/3/s)_. 
Wn + s) = cos(Pis)ui{t n ) H ltj(i n ) 



2Pi(l + tif)J dw 

Here and after we adopt the convention sin(x)/x = 1 when x = 0. 
When n = 0, from the initial condition (12. 2D : i.e., 



— ff(w) sin (A(s - w)) dw. 



5,(0) = («(«)),, 5{(0) = (««),, 
and approximating the integral via the standard trapezoidal rule, we have 



5,(*i) = ooB^rJCtiW), + ^^(««), 



20,(1 + /if) Jo 



T d 



sin (/3,(r - w)) dw 



(2.7), C os( fl ,)^, + ^^,-^^(^)), 



ft 

For n > 1, by the integration by parts, 



2fl(l + Mf) 



thus, 



(2.8) 



— /f H sin (A(s - «;)) dti; = - sin(As)/r(0) + A / /f H cos (A(s - «>)) dw, 



~ i \~ sm (A s )~// 

+ s ) = cos(f3is)ui{t n ) H it,(£ n ) 



2(1 + Mf) 



o 



jf(w) cos (/3,(, - uO) dw - ^^ir(0) 

PI 



Adding (|2.8p with s = r to its evaluation with s = — r, we have the following three-term 

recurrence 

(2.9) 

ui(t n+ i) +ui(t n -i) = 2cos(/?,T)5,(t n )- 1 2 / 

A 1 + Vi) Jo 

In order to design an explicit scheme, we approximate the integral in (|2.9p by the following 
quadrature, 



fPW-fP(-w) cos(A(r-w))dw. 



fPH-fti-w) cos(A(r-w))dw 



(2.10) 



/z»-/r(-r) 



3/r(o)+/r 2 (o)- 4 // n_1 (o) 



TRIGONOMETRIC INTEGRATOR PSEUDOSPECTRAL METHOD FOR SRLW 5 

Here, the first approximation is due to the trapezoidal rule and the second approximation 
is due to a backwards extrapolation. Note that when n = 1, we need an auxiliary approxi- 
mation u(x,t-i) := u(x, — t), which can be achieved similar as (|2.7p from (|2.6p with n = 
and s = — r, 

(2.1D «, - 2^0^), + i^^w,. 

and subjecting (|2.1ip from (|2.7p . we have 

Ui(t-i) + «i(*i) = 2cos(Ar)(n(°)) r 

2.2. Detailed numerical scheme. Let u™ be the approximation of u(xj,t n ) and denote 

by u n the approximation vector with components u". Choosing u° = n^(ajj), the detailed 
numerical scheme computes tJJ as 

(2.12) uj = ]T (^expf^j, j = 0,l,...,M, n = -1, 1, 2, . . . , 

l=-M/2 ^ ' 

where, 

(2.13) (ST), = Oc-^CSO), + (S), - ^^(uOW)),, 

(2.14) (^i), = + 2cos(/3 i r)(S°)) i) 
( n n+i) ; _ -(n^-1)^ + 2cos((3iT)(u n+1 ) l 

(2.15) - , nT 2, W^)l + ~ . n = l,2,..., 

with F n = [(^q) 2 , ("U™) 2 , . . . , (K^-) 2 ] r , and (it n ) ; the discrete Fourier transform coefficients 
for a vector U = [Uq, U±, . . . , Um] T satisfying Uq = Um, defined as 

^ = h E ^ «p ("^) > * = - M A ■ • • » M/2 - 1. 

i=o v 7 

It is readily to see this scheme is fully explicit and easy to implement, with memory cost 
0{M) and computational load 0(M In M) per time step thanks to FFT. 

Remark 2.1. The scheme proposed in this section can be directly extended to the generalized 
symmetric regularized-long-wave equation J2|/, 

utt ~ Uxx + (g{u)) xt ~ u xxtt = 0, 
with g{u) a general nonlinear function. 

3. Numerical results 

In this section numerical results are reported towards testing the accuracy properties 
of the proposed numerical scheme (|2.12p - (|2.15l) . and applying it to study the collisions of 
solitary waves in SRLW. For all the results presented in this section, we always carry out the 
computation on an interval which is large enough such that the periodic boundary conditions 
do not introduce significant aliasing errors relative to the whole space problem. 
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Table 1. Spatial discretization error analysis: discrete Z 2 -error in u at time 
t = 10, 20 under time step r = 0.0001. 


time 


h = 2 h = l h = l/2 


t = 10 
t = 20 


1.5193E-1 1.9505E-4 1.4446E-6 
1.3929E-1 2.0105E-4 3.2758E-6 


Table 2. Temporal discretization error analysis: discrete / 2 -error in u at 
time t = 10, 20 under mesh size h = 1/8. 


time 


r = 0.04 r = 0.02 r = 0.01 r = 0.005 


t = 10 
t = 20 


5.2720E-3 1.3140E-3 3.2810E-4 8.1985E-5 
8.4882E-3 2.0976E-3 5.2155E-4 1.3006E-4 



Table 3. Conservation of invariants analysis: I(u), I(p), E(u, p) and V(u, p) 
at different time levels, with h = 1/2 and r = 0.005. 



time 


J(u) 


«J>) 


E(u,p) 


V(u,p) 


t = 10 


12.0000038 


8.4852811 


27.1529477 


16.8000335 


t = 20 


12.0000081 


8.4852811 


27.1529497 


16.8000348 


t = 30 


12.0000124 


8.4852811 


27.1529515 


16.8000361 


t = 40 


12.0000166 


8.4852811 


27.1529533 


16.8000374 



3.1. Accuracy tests. To test the accuracy properties, we take the initial conditions in 
(|2.2p as the solitary waves (jl.2p with c = y/2, i.e., 



(3.1) u^°\x) = sech 2 ) u^(x) = sech 2 ^"^"^^ tanh • 

To quantify the numerical results, we examine the standard discrete £ 2 -error. 

First, we test the discretization error in space, and in order to do this we take a very fine 
time step r = 0.0001 such that the error from temporal discretization is negligible compared 
to the spatial discretization errors. Tab. [1] lists the errors at time t = 10, 20 for different 
mesh size h. Next, we test the discretization error in time, and again we take a very fine 
mesh size h = 1/8 to eliminate the effect from spatial discretization. Tab. [2] shows the 
errors at time t = 10, 20 for different time step r. Last, we test the conservative properties 
of the invariants given by (jl.5j) - (jl.7p . Tab. [3] lists the computed quantities (ll.5l) - (jl.7p at 
different time levels, under mesh size h = 1/2 and time step r = 0.005. Note that for the 
example presented here, initial p(x,0) is given by u^ ' /y/2 and during computation p(x,t) 
is advanced by the standard midpoint rule with pseudospectral approximations of u x . The 
integrals (ll.5p - (jl.7p are evaluated by the standard trapezoidal quadrature. 

Based on these results, the following are drawn: (i). the scheme (I2.12p -( (2.15p is of spectral- 
order accuracy in space while second-order accuracy in time; (ii). the invariants defined by 
(|1.5p -( fL~7"|) are numerically conserved quite well over a long-time simulation. 



3.2. Applications. Now some examples are reported to study the solitary wave collisions 
in SRLW by applying the proposed numerical scheme. We take the initial conditions in 
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Figure 1. Overtaking collision of two solitary waves with velocity pair {c\ = 
2, c 2 = 6) and x = 12 in d£2J. 
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Figure 2. Contours of overtaking solitary wave collisions. 

(I2T2D as 

(3.2) u (0) (x) = n Cl (x-x ,0) + ii C2 (x + xo,0), u (1) = d t u Cl (x - x , 0) + d t u C2 (x + x , 0), 

and consider various velocities pairs (ci,C2) to cover both cases of overtaking and head-on 
collisions. 

• CASE I: Overtaking solitary wave collisions. 

We consider the initial conditions (13. 2^ with velocities pair {c\,C2) such that 
c\C2 > 0. Fig. [T]depicts the overtaking collision of two solitary waves with velocities 
(ci = 2,C2 = 6) and xq = 12 in (|3.2p . Fig. [2] shows the contours of the evolution 
of overtaking collisions for different initial solitary waves. Based on Figs. [T] and 
[2] and extensive examples, covering a wide range of velocities (01,02), not shown 
here for brevity, overtaking solitary wave collisions are rather elastic with only a 
quite small amount of radiation emitted and no additional visible wave structures 
observed. 

• CASE II: Head-on solitary wave collisions. 

We consider the initial conditions (j3.2j) with velocities pair (01,02) such that 
C1C2 < 0, and study two situations: 

(i) . Initial two solitary waves with equal amplitudes, i.e., c\ = —C2- Fig. [3]depicts 

the head-on collision of two solitary waves with velocities (ci = — 15, C2 = 15) 
and xq = 20 in (|3.2p . It shows that after the collision conspicuous radiations 
are emitted and a number of soliton pairs are produced from either side of 
the radiation. For instance, at time t = 11.8 four soliton pairs are readily 
detectable and at time t = 30.0 five conspicuous soliton pairs are observed. 
Similar results are also observed for other choices of initial equal amplitudes, 
some of which are plotted out in Fig. HI It suggests that there is no upper 
limit on the number of emerging soliton pairs, which indeed appears to be a 
non-decreasing function in the argument of velocity amplitude |ci| (or | C2 1 ) - 

(ii) . Initial two solitary waves with unequal amplitudes, i.e., c\ 7^ — C2. Fig. [5] 

shows the contour and evolution of the head-on collision of two solitary waves 
with velocities {c\ = — 20, C2 = 18) and xq = 20 in (13. 2p . It shows that after 
the collision conspicuous radiations are emitted with an additional soliton pair 
emerging from either side of the radiation. It is interesting to point out that in 
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Figure 3. Head-on collision of two equal-amplitude solitary waves with ve- 
locity pair (ci = —15, C2 = 15) and xq = 20 in (13.2j) . 
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e. = -5, c, = 5, x n = 20 




c. --20, c, = 20,x n = 20 



-100 -50 50 100 




c, =-30, c, = 30, x„ = 20 



-200 -100 100 200 




-200 -100 100 200 



Figure 4. Contours of head-on equal-amplitude solitary wave collisions. 
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Figure 5. Head-on collision of two unequal-amplitude solitary waves with 
velocity pair (ci = —20, ci = 18) and xq = 20 in 
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Figure 6. Head-on collision of two unequal-amplitude solitary waves with 
velocity pair (ci = — 3, C2 = 10) and xq = 20 in (|3.2p . 



this example although the initial amplitudes are nearly equal, the propagation 
result differs much from that of equal-amplitude solitary waves. The head-on 
collision of equal-amplitude solitary waves with (ci = —20, C2 = 20) produces 
four additional soliton pairs after the collision (see Fig. H]), while in Fig. [5] only 
one additional soliton pair emerges. Fig. [U] shows the case (ci = —3, C2 = 10), 
in which the collision is nearly elastic with only a small amount of radiation 
emitted. 

The results presented in this section suggest that the collisions of solitary waves in the 
SRLW equation are rather complicated issues. Also, these results well demonstrate the 
efficiency and high resolution of the proposed scheme (|2.12|) - (j2.15f) for the numerical study 
of the SRLW equation. 



4. Conclusions 

We solved the symmetric regularized-long-wave (SRLW) equation with the utility of pseu- 
dospectral discretization in space combined with a time integration based on trigonomet- 
ric polynomials. The proposed numerical scheme is of second-order accuracy in time and 



12 X. DONG 

spectral-order accuracy in space. The collisions of solitary waves in the SRLW were ex- 
tensively investigated via numerical simulations. The method and results discussed in this 
work are applicable to the physical problem of nonlinear space-charge wave propagation on 
electron beams. 
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